#include <stdio.h>
	// yes, I know stdio.h is not good C++, but I like the *printf()
#include <stdlib.h>
#include <ctype.h>

#define _USE_MATH_DEFINES
#include <math.h>

#ifdef WIN32
#include <windows.h>
#pragma warning(disable:4996)
#endif


// You need to adjust the location of these header files according to your configuration


#include <GL/include/gl.h>
#include <GL/include/glu.h>
#include <GL/include/glut.h>
#include <GL/include/glui.h>

#include "Skeleton.h"

#include <vector>

using namespace std;

//#include "glui.h"

//
//
//	This is a sample OpenGL / GLUT / GLUI program
//
//	The objective is to draw a 3d object and change the color of the axes
//		with radio buttons
//
//	The left mouse button allows rotation
//	The middle mouse button allows scaling
//	The glui window allows:
//		1. The 3d object to be transformed
//		2. The projection to be changed
//		3. The color of the axes to be changed
//		4. The axes to be turned on and off
//		5. The transformations to be reset
//		6. The program to quit
//
//	Author: Joe Graphics
//


//
// constants:
//
// NOTE: There are a bunch of good reasons to use const variables instead
// of #define's.  However, Visual C++ does not allow a const variable
// to be used as an array size or as the case in a switch() statement.  So in
// the following, all constants are const variables except those which need to
// be array sizes or cases in switch() statements.  Those are #defines.
//
//


// This source code has been modified by Guoning Chen since its release


// title of these windows:

const char *WINDOWTITLE = { "OpenGL / GLUT / GLUI Sample -- Joe Graphics" };
const char *GLUITITLE   = { "User Interface Window" };


// what the glui package defines as true and false:

const int GLUITRUE  = { true  };
const int GLUIFALSE = { false };


// the escape key:

#define ESCAPE		0x1b


// initial window size:

const int INIT_WINDOW_SIZE = { 1000 };
const int INIT_WINDOW_WIDTH = { 1000 };
const int INIT_WINDOW_HEIGHT = { 600 };


// size of the box:

const float BOXSIZE = { 2.f };



// multiplication factors for input interaction:
//  (these are known from previous experience)

const float ANGFACT = { 1. };
const float SCLFACT = { 0.005f };


// able to use the left mouse for either rotation or scaling,
// in case have only a 2-button mouse:

enum LeftButton
{
	ROTATE,
	SCALE
};


// minimum allowable scale factor:

const float MINSCALE = { 0.05f };


// active mouse buttons (or them together):

const int LEFT   = { 4 };
const int MIDDLE = { 2 };
const int RIGHT  = { 1 };


// which projection:

enum Projections
{
	ORTHO,
	PERSP
};


// which button:

enum ButtonVals
{
	RESET,
	QUIT,
	LDATA,
	VDATA
};


// window background color (rgba):

const float BACKCOLOR[] = { 0., 0., 0., 0. };


// line width for the axes:

const GLfloat AXES_WIDTH   = { 3. };


// the color numbers:
// this order must match the radio button order

enum Colors
{
	RED,
	YELLOW,
	GREEN,
	CYAN,
	BLUE,
	MAGENTA,
	WHITE,
	BLACK
};


// the color definitions:
// this order must match the radio button order

const GLfloat Colors[8][3] = 
{
	{ 1., 0., 0. },		// red
	{ 1., 1., 0. },		// yellow
	{ 0., 1., 0. },		// green
	{ 0., 1., 1. },		// cyan
	{ 0., 0., 1. },		// blue
	{ 1., 0., 1. },		// magenta
	{ 1., 1., 1. },		// white
	{ 0., 0., 0. },		// black
};


// fog parameters:

const GLfloat FOGCOLOR[4] = { .0, .0, .0, 1. };
const GLenum  FOGMODE     = { GL_LINEAR };
const GLfloat FOGDENSITY  = { 0.30f };
const GLfloat FOGSTART    = { 1.5 };
const GLfloat FOGEND      = { 4. };



//
// non-constant global variables:
//

int	ActiveButton;		// current button that is down
GLuint	AxesList;		// list to hold the axes
int	AxesOn;			// != 0 means to draw the axes
int	DebugOn;			// != 0 means to print debugging info
int	DepthCueOn;		// != 0 means to use intensity depth cueing
GLUI *	glui_subwin;			// instance of glui window
int	GluiWindow;		// the glut id for the glui window
int	LeftButton;		// either ROTATE or SCALE
GLuint	BoxList;		// object display list
int	MainWindow;		// window id for main graphics window
GLfloat	RotMatrix[4][4];	// set by glui rotation widget
float	Scale, Scale2;		// scaling factors
int	WhichColor;		// index into Colors[]
int	WhichProjection;	// ORTHO or PERSP
int	Xmouse, Ymouse;		// mouse values
float	Xrot, Yrot;		// rotation angles in degrees
float	TransXYZ[3];		// set by glui translation widgets


//
// function prototypes:
//

void	Animate( void );
void	Buttons( int );
void	Display( void );
void	DoRasterString( float, float, float, char * );
void	DoStrokeString( float, float, float, float, char * );
float	ElapsedSeconds( void );
void	InitGlui( void );
void	InitGraphics( void );
void	InitLists( void );
void	Keyboard( unsigned char, int, int );
void	MouseButton( int, int, int, int );
void	MouseMotion( int, int );
void	Reset( void );
void	Resize( int, int );
void	Visibility( int );

void	Arrow( float [3], float [3] );
void	Cross( float [3], float [3], float [3] );
float	Dot( float [3], float [3] );
float	Unit( float [3], float [3] );
void	Axes( float );
void	HsvRgb( float[3], float [3] );

void    Display_Model(void);
void    set_view(GLenum mode);
void    set_scene(GLenum mode);
void    display_shape(GLenum mode);
void    Choose_Object();
void	calculateCentroids();

Polyhedron *poly = NULL;

//#dec
//Struct
struct node 
{ 
        float x, y, z;          // location 
		float vx, vy, vz;
		float r, g, b;
		float magnitude;
}; 

struct pointsGrid
{
	vector<vector<node*>> listOfNodes;
};



struct qilVertex
{
	float x,y,z;	//Position
	float nx,ny,nz; //Normal
};
struct qilTriangle{
	int verts[3];
	qilVertex* centroid;
};
struct qilTriangelMesh{
vector<qilVertex*> vertices;
vector<qilTriangle*> triangles;
};

struct vectorsamp{
	qilVertex *tail, *head, *vect;

};
vector<vector<vectorsamp*>> vectorFieldQIL;
qilTriangelMesh *detailedFloor;



//Variables
#define GRIDSIZEX 200
#define GRIDSIZEY 100
#define INTEGRATIONTIME 20

#define AAAA 0.1
#define EPSILON 0.1
#define OMEGA 0.25

#define ALPHA 1.0f
float inst_time = 0.1f;
int inst_timei = 0;
int qilt = 0;
float scaleqil = 2000.;


enum checkBoxValues
{
	VECTORFIELDD,
	VECTORFIELD,
	ADVECTION,
	SPATIALGRADIENTS,
	FTLEFIELD
};
int	VectorFieldOn, AdvectionOn, FTLEFieldOn, SpatialGradientsOn;
int	VectorFieldDOn;
int floorOn;
pointsGrid vectorField;

float flowMapX[INTEGRATIONTIME*10][GRIDSIZEX][GRIDSIZEY];
float flowMapY[INTEGRATIONTIME*10][GRIDSIZEX][GRIDSIZEY];

float SpatailGradientFlowX_dx[INTEGRATIONTIME*10][GRIDSIZEX][GRIDSIZEY];
float SpatailGradientFlowX_dy[INTEGRATIONTIME*10][GRIDSIZEX][GRIDSIZEY];
float SpatailGradientFlowY_dx[INTEGRATIONTIME*10][GRIDSIZEX][GRIDSIZEY];
float SpatailGradientFlowY_dy[INTEGRATIONTIME*10][GRIDSIZEX][GRIDSIZEY];

float FTLEfield[INTEGRATIONTIME*10][GRIDSIZEX][GRIDSIZEY];

unsigned char SpatailGradientFlowX_dxTex[GRIDSIZEX][GRIDSIZEY][4], SpatailGradientFlowX_dyTex[GRIDSIZEX][GRIDSIZEY][4], SpatailGradientFlowY_dxTex[GRIDSIZEX][GRIDSIZEY][4], SpatailGradientFlowY_dyTex[GRIDSIZEX][GRIDSIZEY][4], FTLEFieldTex[GRIDSIZEX][GRIDSIZEY][4];

//functions
void	get_vector_field( float x, float y, float z, float t, float &vxp, float &vyp, float &vzp );
void	get_vector_fieldF( float x, float y, float z, int t, float &vxp, float &vyp, float &vzp ); 
void	create_vector_field(float time);
void	CheckBox( int id );
void	CalculateFlowMaps();
void	CalculateFlowMapsRK4();
void	CalculateSpatialGradients();
void	CalculateFTLEField();
float**	CalculateCauchyGreen(float inputMatrix[][2]);
float	CalcualteEigenVector(float** inputMatrix);
void	genSpatialGradientTexture(int inputTime);	
void	genFTLETexture(int inputTime);	
qilTriangelMesh* loaddata();
void	loadVector();
void	createVectorFieldF(int time);


double radius_factor = 0.9;
int display_mode = 0; 




//extern PlyFile *in_ply;
//extern FILE *this_file;

//
// main program:
//

int
main( int argc, char *argv[] )
{
	// turn on the glut package:
	// (do this before checking argc and argv since it might
	// pull some command line arguments out)

	glutInit( &argc, argv );	


	// setup all the graphics stuff:

	InitGraphics();


	// create the display structures that will not change:

	InitLists();
	create_vector_field(inst_time);
	//CalculateFlowMaps();
	CalculateFlowMapsRK4();
	CalculateSpatialGradients();
	genSpatialGradientTexture(2);	
	CalculateFTLEField();

	// init all the global variables used by Display():
	// this will also post a redisplay
	// it is important to call this before InitGlui()
	// so that the variables that glui will control are correct
	// when each glui widget is created

	Reset();


	// setup all the user interface stuff:

	InitGlui();
	

	// draw the scene once and wait for some interaction:
	// (will never return)

	glutMainLoop();

	// finalize the object if loaded

	

	// this is here to make the compiler happy:

	return 0;
}



//
// this is where one would put code that is to be called
// everytime the glut main loop has nothing to do
//
// this is typically where animation parameters are set
//
// do not call Display() from here -- let glutMainLoop() do it
//

void
Animate( void )
{
	// put animation stuff in here -- change some global variables
	// for Display() to find:
	if(floorOn){

		if(VectorFieldDOn){
			qilt += 1.;
			printf("%d\n",qilt);
			if(qilt > vectorFieldQIL.size()) qilt = 0;
		}

	}
	else{

		if(VectorFieldOn){
			inst_time += 1.;
			create_vector_field(inst_time);
			if(inst_time > INTEGRATIONTIME*10) inst_time = 0.0f;
		}
		// force a call to Display() next time it is convenient:
		if(AdvectionOn){
			inst_timei += 1;	
		
			if(inst_timei >= (INTEGRATIONTIME)*10) inst_timei = 0;
			//if(inst_timei < 10);
			/*getchar()*/
		}
		// force a call to Display() next time it is convenient:
		if(SpatialGradientsOn){
			inst_timei += 1;

			genSpatialGradientTexture((int) inst_timei);
			if(inst_timei >= (INTEGRATIONTIME)*10) inst_timei = 0;
		}
		// force a call to Display() next time it is convenient:
		if(FTLEFieldOn){
			inst_timei += 1;
			genFTLETexture((int) inst_timei);
			if(inst_timei >=  (INTEGRATIONTIME)*10) inst_timei = 0;
		}
	}
	
	glutSetWindow( MainWindow );
	glutPostRedisplay();
}




//
// glui buttons callback:
//

void
Buttons( int id )
{
	switch( id )
	{
		case RESET:
			Reset();
			glui_subwin->sync_live();
			glutSetWindow( MainWindow );
			glutPostRedisplay();
			break;

		case QUIT:
			// gracefully close the glui window:
			// gracefully close out the graphics:
			// gracefully close the graphics window:
			// gracefully exit the program:

			glui_subwin->close();
			glutSetWindow( MainWindow );
			glFinish();
			//glutDestroyWindow( MainWindow );
			exit( 0 );
			break;

		case LDATA:
			detailedFloor = loaddata();
			calculateCentroids();
			floorOn = 1;
			VectorFieldOn = 0;
			AdvectionOn = 0;
			SpatialGradientsOn = 0;
			FTLEFieldOn = 0;
			break;

		case VDATA:
			loadVector();
			createVectorFieldF(10);
			break;

		default:
			fprintf( stderr, "Don't know what to do with Button ID %d\n", id );
	}

}


//
// glui checkbox callback:
//

void
CheckBox( int id )
{
	switch(id){		
	case VECTORFIELD:
		AdvectionOn = 0;
		FTLEFieldOn = 0;
		SpatialGradientsOn = 0;
		break;	

	case VECTORFIELDD:
		VectorFieldDOn = 1;
		AdvectionOn = 0;
		FTLEFieldOn = 0;
		SpatialGradientsOn = 0;
		break;	

	case ADVECTION:
		VectorFieldOn = 0;
		FTLEFieldOn = 0;
		SpatialGradientsOn = 0;
		break;			
		
	case SPATIALGRADIENTS:
		AdvectionOn = 0;
		VectorFieldOn = 0;
		FTLEFieldOn = 0;
		break;	

	case FTLEFIELD:
		VectorFieldOn = 0;
		AdvectionOn = 0;
		SpatialGradientsOn = 0;
		break;		
	}
}


//
// draw the complete scene:
//

void
Display( void )
{
	GLsizei vx, vy, v;		// viewport dimensions
	GLint xl, yb;		// lower-left corner of viewport
	GLfloat scale2;		// real glui scale factor

	if( DebugOn != 0 )
	{
		fprintf( stderr, "Display\n" );
	}


	// set which window we want to do the graphics into:

	glutSetWindow( MainWindow );


	// erase the background:

	glDrawBuffer( GL_BACK );
	glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );
	glEnable( GL_DEPTH_TEST );


	// specify shading to be flat:

	glShadeModel( GL_FLAT );


	// set the viewport to a square centered in the window:

	vx = glutGet( GLUT_WINDOW_WIDTH );
	vy = glutGet( GLUT_WINDOW_HEIGHT );
	v = vx < vy ? vx : vy;			// minimum dimension
	xl = ( vx - v ) / 2;
	yb = ( vy - v ) / 2;
	
	GLUI_Master.get_viewport_area(&xl, &yb, &vx, &vy);
	glViewport( xl, yb,  v, v );
	

	// set the viewing volume:
	// remember that the Z clipping  values are actually
	// given as DISTANCES IN FRONT OF THE EYE
	// USE gluOrtho2D() IF YOU ARE DOING 2D !

	glMatrixMode( GL_PROJECTION );
	glLoadIdentity();
	if( WhichProjection == ORTHO )
		glOrtho( -3., 3.,     -3., 3.,     0.1, 1000. );
	else
		gluPerspective( 90., 1.,	0.1, 1000. );


	// place the objects into the scene:

	glMatrixMode( GL_MODELVIEW );
	glLoadIdentity();


	// set the eye position, look-at position, and up-vector:
	// IF DOING 2D, REMOVE THIS -- OTHERWISE ALL YOUR 2D WILL DISAPPEAR !

	gluLookAt( 0., 0., 3.,     0., 0., 0.,     0., 1., 0. );


	// translate the objects in the scene:
	// note the minus sign on the z value
	// this is to make the appearance of the glui z translate
	// widget more intuitively match the translate behavior
	// DO NOT TRANSLATE IN Z IF YOU ARE DOING 2D !

	glTranslatef( (GLfloat)TransXYZ[0], (GLfloat)TransXYZ[1], -(GLfloat)TransXYZ[2] );


	// rotate the scene:
	// DO NOT ROTATE (EXCEPT ABOUT Z) IF YOU ARE DOING 2D !

	glRotatef( (GLfloat)Yrot, 0., 1., 0. );
	glRotatef( (GLfloat)Xrot, 1., 0., 0. );
	glMultMatrixf( (const GLfloat *) RotMatrix );


	// uniformly scale the scene:

	glScalef( (GLfloat)Scale, (GLfloat)Scale, (GLfloat)Scale );
	scale2 = 1. + Scale2;		// because glui translation starts at 0.
	if( scale2 < MINSCALE )
		scale2 = MINSCALE;
	glScalef( (GLfloat)scale2, (GLfloat)scale2, (GLfloat)scale2 );


	// set the fog parameters:
	// DON'T NEED THIS IF DOING 2D !

	if( DepthCueOn != 0 )
	{
		glFogi( GL_FOG_MODE, FOGMODE );
		glFogfv( GL_FOG_COLOR, FOGCOLOR );
		glFogf( GL_FOG_DENSITY, FOGDENSITY );
		glFogf( GL_FOG_START, FOGSTART );
		glFogf( GL_FOG_END, FOGEND );
		glEnable( GL_FOG );
	}
	else
	{
		glDisable( GL_FOG );
	}

	// Let us disable lighting right now
	glDisable(GL_LIGHTING);


	// possibly draw the axes:

	if( AxesOn != 0 )
	{
		glColor3fv( &Colors[WhichColor][0] );
		glCallList( AxesList );
	}


	// set the color of the object:

	glColor3fv( Colors[WhichColor] );


	// draw the current object:

	//glCallList( BoxList );


	// draw some gratuitous text that just rotates on top of the scene:

	glDisable( GL_DEPTH_TEST );
	glColor3f( 0., 1., 1. );
	//DoRasterString( 0., 1., 0., "Text That Moves" );


	// draw some gratuitous text that is fixed on the screen:
	// the projection matrix is reset to define a scene whose
	// world coordinate system goes from 0-100 in each axis
	// this is called "percent units", and is just a convenience
	// the modelview matrix is reset to identity as we don't
	// want to transform these coordinates

	glDisable( GL_DEPTH_TEST );
	glMatrixMode( GL_PROJECTION );
	glLoadIdentity();
	gluOrtho2D( 0., 100.,     0., 100. );
	glMatrixMode( GL_MODELVIEW );
	glLoadIdentity();
	glColor3f( 1., 1., 1. );
	//DoRasterString( 5., 5., 0., "t = " + itoa(inst_time) );
	
	char str[20];
	if(VectorFieldOn) sprintf(str, "%f", (float)inst_time/10.);
	else sprintf(str, "%f", (float)inst_timei/10.);
	DoRasterString( 5., 5., 0., str );

	// Render the loaded object
	set_view(GL_RENDER);
	
	glTranslatef(0.0, 0.0, -3.0);

	glTranslatef( (GLfloat)TransXYZ[0], (GLfloat)TransXYZ[1], -(GLfloat)TransXYZ[2] );
	
	

	glScalef( (GLfloat)scale2, (GLfloat)scale2, (GLfloat)scale2 );
	
	glRotatef( (GLfloat)Yrot, 0., 1., 0. );
	glRotatef( (GLfloat)Xrot, 1., 0., 0. );
	glMultMatrixf( (const GLfloat *) RotMatrix );
	

	

	display_shape(GL_RENDER);

	Sleep(250);

	// swap the double-buffered framebuffers:

	glutSwapBuffers();


	// be sure the graphics buffer has been sent:
	// note: be sure to use glFlush() here, not glFinish() !

	glFlush();
}



//
// use glut to display a string of characters using a raster font:
//

void
DoRasterString( float x, float y, float z, char *s )
{
	char c;			// one character to print

	glRasterPos3f( (GLfloat)x, (GLfloat)y, (GLfloat)z );
	for( ; ( c = *s ) != '\0'; s++ )
	{
		glutBitmapCharacter( GLUT_BITMAP_TIMES_ROMAN_24, c );
	}
}



//
// use glut to display a string of characters using a stroke font:
//

void
DoStrokeString( float x, float y, float z, float ht, char *s )
{
	char c;			// one character to print
	float sf;		// the scale factor

	glPushMatrix();
		glTranslatef( (GLfloat)x, (GLfloat)y, (GLfloat)z );
		sf = ht / ( 119.05 + 33.33 );
		glScalef( (GLfloat)sf, (GLfloat)sf, (GLfloat)sf );
		for( ; ( c = *s ) != '\0'; s++ )
		{
			glutStrokeCharacter( GLUT_STROKE_ROMAN, c );
		}
	glPopMatrix();
}



//
// return the number of seconds since the start of the program:
//

float
ElapsedSeconds( void )
{
	// get # of milliseconds since the start of the program:

	int ms = glutGet( GLUT_ELAPSED_TIME );

	// convert it to seconds:

	return (float)ms / 1000.;
}



//
// initialize the glui window:
//#IGF

void
InitGlui( void )
{
	GLUI_Panel *panel;
	GLUI_RadioGroup *group;
	GLUI_Rotation *rot;
	GLUI_Translation *trans, *scale;


	// setup the glui window:

	/*glutInitWindowPosition( INIT_WINDOW_WIDTH + 50, 0 );
	Glui = GLUI_Master.create_glui( (char *) GLUITITLE );*/

	//Adding the graphics display window to the glui window
	GLUI_Master.set_glutKeyboardFunc( Keyboard );
	GLUI_Master.set_glutMouseFunc ( MouseButton );
	GLUI_Master.set_glutReshapeFunc( Resize );

	glui_subwin =  GLUI_Master.create_glui_subwindow(MainWindow, GLUI_SUBWINDOW_RIGHT);
	glui_subwin->set_main_gfx_window( MainWindow );

	glui_subwin->add_statictext( (char *) GLUITITLE );
	glui_subwin->add_separator();

	glui_subwin->add_checkbox( "Axes", &AxesOn );

	glui_subwin->add_checkbox( "Perspective", &WhichProjection );

	glui_subwin->add_checkbox( "Intensity Depth Cue", &DepthCueOn );

	
	panel = glui_subwin->add_panel( "Object Transformation" );

		rot = glui_subwin->add_rotation_to_panel( panel, "Rotation", (float *) RotMatrix );

		// allow the object to be spun via the glui rotation widget:

		rot->set_spin( 1.0 );


		glui_subwin->add_column_to_panel( panel, GLUIFALSE );
		scale = glui_subwin->add_translation_to_panel( panel, "Scale",  GLUI_TRANSLATION_Y , &Scale2 );
		scale->set_speed( 0.005f );

		glui_subwin->add_column_to_panel( panel, GLUIFALSE );
		trans = glui_subwin->add_translation_to_panel( panel, "Trans XY", GLUI_TRANSLATION_XY, &TransXYZ[0] );
		trans->set_speed( 0.05f );

		glui_subwin->add_column_to_panel( panel, GLUIFALSE );
		trans = glui_subwin->add_translation_to_panel( panel, "Trans Z",  GLUI_TRANSLATION_Z , &TransXYZ[2] );
		trans->set_speed( 0.05f );

		glui_subwin->add_checkbox("Vector", &VectorFieldOn, VECTORFIELD, (GLUI_Update_CB) CheckBox);
		glui_subwin->add_checkbox("Advection", &AdvectionOn, ADVECTION, (GLUI_Update_CB) CheckBox);
		glui_subwin->add_checkbox("Spatial Gradients", &SpatialGradientsOn, SPATIALGRADIENTS, (GLUI_Update_CB) CheckBox);
		glui_subwin->add_checkbox("FTLE", &FTLEFieldOn, FTLEFIELD, (GLUI_Update_CB) CheckBox);

	glui_subwin->add_button("Load Floor", LDATA, (GLUI_Update_CB) Buttons);
	glui_subwin->add_button("Load Vector", VDATA, (GLUI_Update_CB) Buttons);

	glui_subwin->add_checkbox("Vector D", &VectorFieldDOn, VECTORFIELDD, (GLUI_Update_CB) CheckBox);

	glui_subwin->add_checkbox( "Debug", &DebugOn );


	panel = glui_subwin->add_panel( "", GLUIFALSE );

	glui_subwin->add_button_to_panel( panel, "Reset", RESET, (GLUI_Update_CB) Buttons );

	glui_subwin->add_column_to_panel( panel, GLUIFALSE );

	glui_subwin->add_button_to_panel( panel, "Quit", QUIT, (GLUI_Update_CB) Buttons );


	// tell glui what graphics window it needs to post a redisplay to:

	glui_subwin->set_main_gfx_window( MainWindow );


	// set the graphics window's idle function:

	//GLUI_Master.set_glutIdleFunc( NULL );
	GLUI_Master.set_glutIdleFunc( Animate );
}



//
// initialize the glut and OpenGL libraries:
//	also setup display lists and callback functions
//

void
InitGraphics( void )
{
	// setup the display mode:
	// ( *must* be done before call to glutCreateWindow() )
	// ask for color, double-buffering, and z-buffering:

	glutInitDisplayMode( GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH );


	// set the initial window configuration:

	glutInitWindowPosition( 0, 0 );
	glutInitWindowSize( INIT_WINDOW_WIDTH, INIT_WINDOW_HEIGHT );


	// open the window and set its title:

	MainWindow = glutCreateWindow( WINDOWTITLE );
	glutSetWindowTitle( WINDOWTITLE );


	// setup the clear values:

	glClearColor( BACKCOLOR[0], BACKCOLOR[1], BACKCOLOR[2], BACKCOLOR[3] );


	// setup the callback routines:


	// DisplayFunc -- redraw the window
	// ReshapeFunc -- handle the user resizing the window
	// KeyboardFunc -- handle a keyboard input
	// MouseFunc -- handle the mouse button going down or up
	// MotionFunc -- handle the mouse moving with a button down
	// PassiveMotionFunc -- handle the mouse moving with a button up
	// VisibilityFunc -- handle a change in window visibility
	// EntryFunc	-- handle the cursor entering or leaving the window
	// SpecialFunc -- handle special keys on the keyboard
	// SpaceballMotionFunc -- handle spaceball translation
	// SpaceballRotateFunc -- handle spaceball rotation
	// SpaceballButtonFunc -- handle spaceball button hits
	// ButtonBoxFunc -- handle button box hits
	// DialsFunc -- handle dial rotations
	// TabletMotionFunc -- handle digitizing tablet motion
	// TabletButtonFunc -- handle digitizing tablet button hits
	// MenuStateFunc -- declare when a pop-up menu is in use
	// TimerFunc -- trigger something to happen a certain time from now
	// IdleFunc -- what to do when nothing else is going on

	glutSetWindow( MainWindow );
	glutDisplayFunc( Display );
	//glutDisplayFunc( Display_Model );
	glutReshapeFunc( Resize );
	glutKeyboardFunc( Keyboard );
	glutMouseFunc( MouseButton );
	glutMotionFunc( MouseMotion );
	glutPassiveMotionFunc( NULL );
	glutVisibilityFunc( Visibility );
	glutEntryFunc( NULL );
	glutSpecialFunc( NULL );
	glutSpaceballMotionFunc( NULL );
	glutSpaceballRotateFunc( NULL );
	glutSpaceballButtonFunc( NULL );
	glutButtonBoxFunc( NULL );
	glutDialsFunc( NULL );
	glutTabletMotionFunc( NULL );
	glutTabletButtonFunc( NULL );
	glutMenuStateFunc( NULL );
	glutTimerFunc( 0, NULL, 0 );

	// DO NOT SET THE GLUT IDLE FUNCTION HERE !!
	// glutIdleFunc( NULL );
	// let glui take care of it in InitGlui()
}




//
// initialize the display lists that will not change:
//

void
InitLists( void )
{
	float dx = BOXSIZE / 2.;
	float dy = BOXSIZE / 2.;
	float dz = BOXSIZE / 2.;

	// create the object:

	BoxList = glGenLists( 1 );
	glNewList( BoxList, GL_COMPILE );

		glBegin( GL_QUADS );

			glColor3f( 0., 0., 1. );
			glNormal3f( 0., 0.,  1. );
				glVertex3f( -dx, -dy,  dz );
				glVertex3f(  dx, -dy,  dz );
				glVertex3f(  dx,  dy,  dz );
				glVertex3f( -dx,  dy,  dz );

			glNormal3f( 0., 0., -1. );
				glTexCoord2f( 0., 0. );
				glVertex3f( -dx, -dy, -dz );
				glTexCoord2f( 0., 1. );
				glVertex3f( -dx,  dy, -dz );
				glTexCoord2f( 1., 1. );
				glVertex3f(  dx,  dy, -dz );
				glTexCoord2f( 1., 0. );
				glVertex3f(  dx, -dy, -dz );

			glColor3f( 1., 0., 0. );
			glNormal3f(  1., 0., 0. );
				glVertex3f(  dx, -dy,  dz );
				glVertex3f(  dx, -dy, -dz );
				glVertex3f(  dx,  dy, -dz );
				glVertex3f(  dx,  dy,  dz );

			glNormal3f( -1., 0., 0. );
				glVertex3f( -dx, -dy,  dz );
				glVertex3f( -dx,  dy,  dz );
				glVertex3f( -dx,  dy, -dz );
				glVertex3f( -dx, -dy, -dz );

			glColor3f( 0., 1., 0. );
			glNormal3f( 0.,  1., 0. );
				glVertex3f( -dx,  dy,  dz );
				glVertex3f(  dx,  dy,  dz );
				glVertex3f(  dx,  dy, -dz );
				glVertex3f( -dx,  dy, -dz );

			glNormal3f( 0., -1., 0. );
				glVertex3f( -dx, -dy,  dz );
				glVertex3f( -dx, -dy, -dz );
				glVertex3f(  dx, -dy, -dz );
				glVertex3f(  dx, -dy,  dz );

		glEnd();

	glEndList();


	// create the axes:

	AxesList = glGenLists( 1 );
	glNewList( AxesList, GL_COMPILE );
		glLineWidth( AXES_WIDTH );
			Axes( 1.5 );
		glLineWidth( 1. );
	glEndList();
}



//
// the keyboard callback:
//

void
Keyboard( unsigned char c, int x, int y )
{
	if( DebugOn != 0 )
		fprintf( stderr, "Keyboard: '%c' (0x%0x)\n", c, c );

	switch( c )
	{
		case 'o':
		case 'O':
			WhichProjection = ORTHO;
			break;

		case 'p':
		case 'P':
			WhichProjection = PERSP;
			break;

		case 'q':
		case 'Q':
		case ESCAPE:
			Buttons( QUIT );	// will not return here
			break;			// happy compiler

		case 'r':
		case 'R':
			LeftButton = ROTATE;
			break;

		case 's':
		case 'S':
			LeftButton = SCALE;
			break;

		default:
			fprintf( stderr, "Don't know what to do with keyboard hit: '%c' (0x%0x)\n", c, c );
	}


	// synchronize the GLUI display with the variables:

	glui_subwin->sync_live();


	// force a call to Display():

	glutSetWindow( MainWindow );
	glutPostRedisplay();
}



//
// called when the mouse button transitions down or up:
//

void
MouseButton( int button, int state, int x, int y )
{
	int b;			// LEFT, MIDDLE, or RIGHT

	if( DebugOn != 0 )
		fprintf( stderr, "MouseButton: %d, %d, %d, %d\n", button, state, x, y );

	
	// get the proper button bit mask:

	switch( button )
	{
		case GLUT_LEFT_BUTTON:
			b = LEFT;		break;

		case GLUT_MIDDLE_BUTTON:
			b = MIDDLE;		break;

		case GLUT_RIGHT_BUTTON:
			b = RIGHT;		break;

		default:
			b = 0;
			fprintf( stderr, "Unknown mouse button: %d\n", button );
	}


	// button down sets the bit, up clears the bit:

	if( state == GLUT_DOWN )
	{
		Xmouse = x;
		Ymouse = y;
		ActiveButton |= b;		// set the proper bit
	}
	else
	{
		ActiveButton &= ~b;		// clear the proper bit
	}
}



//
// called when the mouse moves while a button is down:
//

void
MouseMotion( int x, int y )
{
	int dx, dy;		// change in mouse coordinates

	if( DebugOn != 0 )
		fprintf( stderr, "MouseMotion: %d, %d\n", x, y );


	dx = x - Xmouse;		// change in mouse coords
	dy = y - Ymouse;

	if( ( ActiveButton & LEFT ) != 0 )
	{
		switch( LeftButton )
		{
			case ROTATE:
				Xrot += ( ANGFACT*dy );
				Yrot += ( ANGFACT*dx );
				break;

			case SCALE:
				Scale += SCLFACT * (float) ( dx - dy );
				if( Scale < MINSCALE )
					Scale = MINSCALE;
				break;
		}
	}


	if( ( ActiveButton & MIDDLE ) != 0 )
	{
		Scale += SCLFACT * (float) ( dx - dy );

		// keep object from turning inside-out or disappearing:

		if( Scale < MINSCALE )
			Scale = MINSCALE;
	}

	Xmouse = x;			// new current position
	Ymouse = y;

	glutSetWindow( MainWindow );
	glutPostRedisplay();
}



//
// reset the transformations and the colors:
//
// this only sets the global variables --
// the glut main loop is responsible for redrawing the scene
//

void
Reset( void )
{
	ActiveButton = 0;
	AxesOn = GLUITRUE;
	DebugOn = GLUIFALSE;
	DepthCueOn = GLUIFALSE;
	LeftButton = ROTATE;
	Scale  = 1.0;
	Scale2 = 0.0;		// because we add 1. to it in Display()
	WhichColor = WHITE;
	WhichProjection = PERSP;
	Xrot = Yrot = 0.;
	TransXYZ[0] = TransXYZ[1] = TransXYZ[2] = 0.;

	                  RotMatrix[0][1] = RotMatrix[0][2] = RotMatrix[0][3] = 0.;
	RotMatrix[1][0]                   = RotMatrix[1][2] = RotMatrix[1][3] = 0.;
	RotMatrix[2][0] = RotMatrix[2][1]                   = RotMatrix[2][3] = 0.;
	RotMatrix[3][0] = RotMatrix[3][1] = RotMatrix[3][3]                   = 0.;
	RotMatrix[0][0] = RotMatrix[1][1] = RotMatrix[2][2] = RotMatrix[3][3] = 1.;

	inst_time = 0.; inst_timei = 0;
}



//
// called when user resizes the window:
//

void
Resize( int width, int height )
{
	if( DebugOn != 0 )
		fprintf( stderr, "ReSize: %d, %d\n", width, height );

	// don't really need to do anything since window size is
	// checked each time in Display():

	glutSetWindow( MainWindow );
	glutPostRedisplay();
}


//
// handle a change to the window's visibility:
//

void
Visibility ( int state )
{
	if( DebugOn != 0 )
		fprintf( stderr, "Visibility: %d\n", state );

	if( state == GLUT_VISIBLE )
	{
		glutSetWindow( MainWindow );
		glutPostRedisplay();
	}
	else
	{
		// could optimize by keeping track of the fact
		// that the window is not visible and avoid
		// animating or redrawing it ...
	}
}




//////////////////////////////////////////  EXTRA HANDY UTILITIES:  /////////////////////////////

// size of wings as fraction of length:

#define WINGS	0.10


// axes:

#define X	1
#define Y	2
#define Z	3


// x, y, z, axes:

static float axx[3] = { 1., 0., 0. };
static float ayy[3] = { 0., 1., 0. };
static float azz[3] = { 0., 0., 1. };


void
Arrow( float tail[3], float head[3] )
{
	float u[3], v[3], w[3];		// arrow coordinate system
	float d;			// wing distance
	float x, y, z;			// point to plot
	float mag;			// magnitude of major direction
	float f;			// fabs of magnitude
	int axis;			// which axis is the major


	// set w direction in u-v-w coordinate system:

	w[0] = head[0] - tail[0];
	w[1] = head[1] - tail[1];
	w[2] = head[2] - tail[2];


	// determine major direction:

	axis = X;
	mag = fabs( w[0] );
	if( (f=fabs(w[1]))  > mag )
	{
		axis = Y;
		mag = f;
	}
	if( (f=fabs(w[2]))  > mag )
	{
		axis = Z;
		mag = f;
	}


	// set size of wings and turn w into a Unit vector:

	d = WINGS * Unit( w, w );


	// draw the shaft of the arrow:

	glBegin( GL_LINE_STRIP );
		glVertex3fv( tail );
		glVertex3fv( head );
	glEnd();

	// draw two sets of wings in the non-major directions:

	if( axis != X )
	{
		Cross( w, axx, v );
		(void) Unit( v, v );
		Cross( v, w, u  );
		x = head[0] + d * ( u[0] - w[0] );
		y = head[1] + d * ( u[1] - w[1] );
		z = head[2] + d * ( u[2] - w[2] );
		glBegin( GL_LINE_STRIP );
			glVertex3fv( head );
			glVertex3f( x, y, z );
		glEnd();
		x = head[0] + d * ( -u[0] - w[0] );
		y = head[1] + d * ( -u[1] - w[1] );
		z = head[2] + d * ( -u[2] - w[2] );
		glBegin( GL_LINE_STRIP );
			glVertex3fv( head );
			glVertex3f( x, y, z );
		glEnd();
	}


	if( axis != Y )
	{
		Cross( w, ayy, v );
		(void) Unit( v, v );
		Cross( v, w, u  );
		x = head[0] + d * ( u[0] - w[0] );
		y = head[1] + d * ( u[1] - w[1] );
		z = head[2] + d * ( u[2] - w[2] );
		glBegin( GL_LINE_STRIP );
			glVertex3fv( head );
			glVertex3f( x, y, z );
		glEnd();
		x = head[0] + d * ( -u[0] - w[0] );
		y = head[1] + d * ( -u[1] - w[1] );
		z = head[2] + d * ( -u[2] - w[2] );
		glBegin( GL_LINE_STRIP );
			glVertex3fv( head );
			glVertex3f( x, y, z );
		glEnd();
	}



	if( axis != Z )
	{
		Cross( w, azz, v );
		(void) Unit( v, v );
		Cross( v, w, u  );
		x = head[0] + d * ( u[0] - w[0] );
		y = head[1] + d * ( u[1] - w[1] );
		z = head[2] + d * ( u[2] - w[2] );
		glBegin( GL_LINE_STRIP );
			glVertex3fv( head );
			glVertex3f( x, y, z );
		glEnd();
		x = head[0] + d * ( -u[0] - w[0] );
		y = head[1] + d * ( -u[1] - w[1] );
		z = head[2] + d * ( -u[2] - w[2] );
		glBegin( GL_LINE_STRIP );
			glVertex3fv( head );
			glVertex3f( x, y, z );
		glEnd();
	}
}



float
Dot( float v1[3], float v2[3] )
{
	return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
}



void
Cross( float v1[3], float v2[3], float vout[3] )
{
	float tmp[3];

	tmp[0] = v1[1]*v2[2] - v2[1]*v1[2];
	tmp[1] = v2[0]*v1[2] - v1[0]*v2[2];
	tmp[2] = v1[0]*v2[1] - v2[0]*v1[1];

	vout[0] = tmp[0];
	vout[1] = tmp[1];
	vout[2] = tmp[2];
}



float
Unit( float vin[3], float vout[3] )
{
	float dist, f ;

	dist = vin[0]*vin[0] + vin[1]*vin[1] + vin[2]*vin[2];

	if( dist > 0.0 )
	{
		dist = sqrt( dist );
		f = 1. / dist;
		vout[0] = f * vin[0];
		vout[1] = f * vin[1];
		vout[2] = f * vin[2];
	}
	else
	{
		vout[0] = vin[0];
		vout[1] = vin[1];
		vout[2] = vin[2];
	}

	return dist;
}



// the stroke characters 'X' 'Y' 'Z' :

static float xx[] = {
		0.f, 1.f, 0.f, 1.f
	      };

static float xy[] = {
		-.5f, .5f, .5f, -.5f
	      };

static int xorder[] = {
		1, 2, -3, 4
		};


static float yx[] = {
		0.f, 0.f, -.5f, .5f
	      };

static float yy[] = {
		0.f, .6f, 1.f, 1.f
	      };

static int yorder[] = {
		1, 2, 3, -2, 4
		};


static float zx[] = {
		1.f, 0.f, 1.f, 0.f, .25f, .75f
	      };

static float zy[] = {
		.5f, .5f, -.5f, -.5f, 0.f, 0.f
	      };

static int zorder[] = {
		1, 2, 3, 4, -5, 6
		};


// fraction of the length to use as height of the characters:

const float LENFRAC = 0.10f;


// fraction of length to use as start location of the characters:

const float BASEFRAC = 1.10f;


//
//	Draw a set of 3D axes:
//	(length is the axis length in world coordinates)
//

void
Axes( float length )
{
	int i, j;			// counters
	float fact;			// character scale factor
	float base;			// character start location


	glBegin( GL_LINE_STRIP );
		glVertex3f( length, 0., 0. );
		glVertex3f( 0., 0., 0. );
		glVertex3f( 0., length, 0. );
	glEnd();
	glBegin( GL_LINE_STRIP );
		glVertex3f( 0., 0., 0. );
		glVertex3f( 0., 0., length );
	glEnd();

	fact = LENFRAC * length;
	base = BASEFRAC * length;

	glBegin( GL_LINE_STRIP );
		for( i = 0; i < 4; i++ )
		{
			j = xorder[i];
			if( j < 0 )
			{
				
				glEnd();
				glBegin( GL_LINE_STRIP );
				j = -j;
			}
			j--;
			glVertex3f( base + fact*xx[j], fact*xy[j], 0.0 );
		}
	glEnd();

	glBegin( GL_LINE_STRIP );
		for( i = 0; i < 5; i++ )
		{
			j = yorder[i];
			if( j < 0 )
			{
				
				glEnd();
				glBegin( GL_LINE_STRIP );
				j = -j;
			}
			j--;
			glVertex3f( fact*yx[j], base + fact*yy[j], 0.0 );
		}
	glEnd();

	glBegin( GL_LINE_STRIP );
		for( i = 0; i < 6; i++ )
		{
			j = zorder[i];
			if( j < 0 )
			{
				
				glEnd();
				glBegin( GL_LINE_STRIP );
				j = -j;
			}
			j--;
			glVertex3f( 0.0, fact*zy[j], base + fact*zx[j] );
		}
	glEnd();

}




//
// routine to convert HSV to RGB
//
// Reference:  Foley, van Dam, Feiner, Hughes,
//		"Computer Graphics Principles and Practices,"
//		Additon-Wesley, 1990, pp592-593.


void
HsvRgb( float hsv[3], float rgb[3] )
{
	float h, s, v;			// hue, sat, value
	float r, g, b;			// red, green, blue
	float i, f, p, q, t;		// interim values


	// guarantee valid input:

	h = hsv[0] / 60.;
	while( h >= 6. )	h -= 6.;
	while( h <  0. ) 	h += 6.;

	s = hsv[1];
	if( s < 0. )
		s = 0.;
	if( s > 1. )
		s = 1.;

	v = hsv[2];
	if( v < 0. )
		v = 0.;
	if( v > 1. )
		v = 1.;


	// if sat==0, then is a gray:

	if( s == 0.0 )
	{
		rgb[0] = rgb[1] = rgb[2] = v;
		return;
	}


	// get an rgb from the hue itself:
	
	i = floor( h );
	f = h - i;
	p = v * ( 1. - s );
	q = v * ( 1. - s*f );
	t = v * ( 1. - ( s * (1.-f) ) );

	switch( (int) i )
	{
		case 0:
			r = v;	g = t;	b = p;
			break;
	
		case 1:
			r = q;	g = v;	b = p;
			break;
	
		case 2:
			r = p;	g = v;	b = t;
			break;
	
		case 3:
			r = p;	g = q;	b = v;
			break;
	
		case 4:
			r = t;	g = p;	b = v;
			break;
	
		case 5:
			r = v;	g = p;	b = q;
			break;
	}


	rgb[0] = r;
	rgb[1] = g;
	rgb[2] = b;
}


void set_view(GLenum mode)
{
	icVector3 up, ray, view;
	GLfloat light_ambient0[] = { 0.3, 0.3, 0.3, 1.0 };
	GLfloat light_diffuse0[] = { 0.7, 0.7, 0.7, 1.0 };
	GLfloat light_specular0[] = { 0.0, 0.0, 0.0, 1.0 };
	GLfloat light_ambient1[] = { 0.0, 0.0, 0.0, 1.0 };
	GLfloat light_diffuse1[] = { 0.5, 0.5, 0.5, 1.0 };
	GLfloat light_specular1[] = { 0.0, 0.0, 0.0, 1.0 };
	GLfloat light_ambient2[] = { 1.0, 1.0, 1.0, 1.0 };
	GLfloat light_diffuse2[] = { 1.0, 1.0, 1.0, 1.0 };
	GLfloat light_specular2[] = { 1.0, 1.0, 1.0, 1.0 };
	GLfloat light_position[] = { 0.0, 0.0, 0.0, 1.0 };

	glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient0);
	glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse0);
	glLightfv(GL_LIGHT0, GL_SPECULAR, light_specular0);
	glLightfv(GL_LIGHT1, GL_AMBIENT, light_ambient1);
	glLightfv(GL_LIGHT1, GL_DIFFUSE, light_diffuse1);
	glLightfv(GL_LIGHT1, GL_SPECULAR, light_specular1);


  glMatrixMode(GL_PROJECTION);
	if (mode == GL_RENDER)
		glLoadIdentity();

	if (WhichProjection == ORTHO)
		glOrtho(-radius_factor, radius_factor, -radius_factor, radius_factor, 0.0, 40.0);
	else
		gluPerspective(45.0, 1.0, 0.1, 40.0);

	glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
	light_position[0] = 5.5;
	light_position[1] = 0.0;
	light_position[2] = 0.0;
	glLightfv(GL_LIGHT0, GL_POSITION, light_position);
	light_position[0] = -0.1;
	light_position[1] = 0.0;
	light_position[2] = 0.0;
	glLightfv(GL_LIGHT2, GL_POSITION, light_position);
}

void set_scene(GLenum mode)
{
	
	
}


void display_shape(GLenum mode)
{

	if(floorOn){
		glColor3f(1., 0., 0.);					
		glBegin(GL_LINES);
		for (int i = 0; i< detailedFloor->triangles.size() ; i++){
			for (int j = 0; j< detailedFloor->vertices.size() ; j++){				
				glVertex3f(detailedFloor->vertices[detailedFloor->triangles[i]->verts[0]]->x/scaleqil, 0., detailedFloor->vertices[detailedFloor->triangles[i]->verts[0]]->z/scaleqil);
				glVertex3f(detailedFloor->vertices[detailedFloor->triangles[i]->verts[1]]->x/scaleqil, 0., detailedFloor->vertices[detailedFloor->triangles[i]->verts[1]]->z/scaleqil);
				
				glVertex3f(detailedFloor->vertices[detailedFloor->triangles[i]->verts[1]]->x/scaleqil, 0., detailedFloor->vertices[detailedFloor->triangles[i]->verts[1]]->z/scaleqil);
				glVertex3f(detailedFloor->vertices[detailedFloor->triangles[i]->verts[2]]->x/scaleqil, 0., detailedFloor->vertices[detailedFloor->triangles[i]->verts[2]]->z/scaleqil);
								
				glVertex3f(detailedFloor->vertices[detailedFloor->triangles[i]->verts[2]]->x/scaleqil, 0., detailedFloor->vertices[detailedFloor->triangles[i]->verts[2]]->z/scaleqil);
				glVertex3f(detailedFloor->vertices[detailedFloor->triangles[i]->verts[0]]->x/scaleqil, 0., detailedFloor->vertices[detailedFloor->triangles[i]->verts[0]]->z/scaleqil);
						
			}
		}
		glEnd();
		if(VectorFieldDOn){
			//glBegin(GL_LINES);	
			///*for (int tim = 0; tim < 40; tim++){qilt = tim;*/
			//		
			//	for (int i = 0; i< vectorFieldQIL[qilt].size(); i++){
			//		
			//		glColor3f(0., 1., 0.);	
			//		printf("%f, %f, %f\n", vectorFieldQIL[qilt][i]->tail->x/scaleqil, vectorFieldQIL[qilt][i]->tail->y/scaleqil, vectorFieldQIL[qilt][i]->tail->z/scaleqil);
			//		float tail[3] = {vectorFieldQIL[qilt][i]->tail->x/scaleqil, vectorFieldQIL[qilt][i]->tail->y/scaleqil, vectorFieldQIL[qilt][i]->tail->z/scaleqil};					
			//		float head[3] = {vectorFieldQIL[qilt][i]->head->x/scaleqil, vectorFieldQIL[qilt][i]->head->y/scaleqil, vectorFieldQIL[qilt][i]->head->z/scaleqil};
			//		Arrow( tail,head );				
			//	
			//	}
			//glEnd();	
			////}

			glBegin(GL_LINES);	
			for (int i = 0; i< GRIDSIZEY; i++){
				for (int j = 0; j< GRIDSIZEX; j++){				
						/*glColor3f(vectorField.listOfNodes[i][j]->r, vectorField.listOfNodes[i][j]->g, vectorField.listOfNodes[i][j]->b);					*/
					if(vectorField.listOfNodes[i][j]->magnitude > 0.){
						glColor3f(0., 1., 0.);	
						float tail[3] = {vectorField.listOfNodes[i][j]->x, 0., vectorField.listOfNodes[i][j]->z};					
						float head[3] = {vectorField.listOfNodes[i][j]->x + vectorField.listOfNodes[i][j]->vx, 0. , vectorField.listOfNodes[i][j]->z + vectorField.listOfNodes[i][j]->vz};
						Arrow( tail,head);				
					}
				}
			}
			glEnd();	
		}
	}
	else{
		unsigned int i, j;
		GLfloat mat_diffuse[4];
		int filter;
		//glColor3f(0.0f, 0.0f, 0.0f);
		if(VectorFieldOn){
			glBegin(GL_LINES);	
			for (int i = 0; i< GRIDSIZEX; i++){
				for (int j = 0; j< GRIDSIZEY; j++){				
						glColor3f(vectorField.listOfNodes[i][j]->r, vectorField.listOfNodes[i][j]->g, vectorField.listOfNodes[i][j]->b);					
						float tail[3] = {vectorField.listOfNodes[i][j]->x, vectorField.listOfNodes[i][j]->y, vectorField.listOfNodes[i][j]->z};					
						float head[3] = {vectorField.listOfNodes[i][j]->x + (vectorField.listOfNodes[i][j]->vx/(/*vectorField.listOfNodes[i][j]->magnitude**/GRIDSIZEX)), vectorField.listOfNodes[i][j]->y + (vectorField.listOfNodes[i][j]->vy/(/*vectorField.listOfNodes[i][j]->magnitude**/GRIDSIZEY)), vectorField.listOfNodes[i][j]->z + (vectorField.listOfNodes[i][j]->vz/(/*vectorField.listOfNodes[i][j]->magnitude**/GRIDSIZEX))};
						Arrow( tail,head);				
				}
			}
			glEnd();	
		}
		//glColor3f(0.0f, 0.0f, 0.0f);
		if(AdvectionOn){
		
			glBegin(GL_POINTS);	
			for (int i = 0; i< GRIDSIZEX; i++){
				for (int j = 0; j< GRIDSIZEY; j++){				
					glColor3f((float)i/(float)GRIDSIZEX, (float)j/(float)GRIDSIZEY, 0.);					
					glVertex3f(flowMapX[inst_timei][i][j], flowMapY[inst_timei][i][j], 0.);
					
				}
			}
			glEnd();	
		}
		//Display spatial gradients as textures
		if(SpatialGradientsOn){
			glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP ); 
			glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP ); 
			glTexEnvf( GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE ); 
			/*filter = GL_NEAREST; */
			filter = GL_LINEAR; 
		

			glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, filter ); 
			glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, filter ); 

			glPixelStorei( GL_UNPACK_ALIGNMENT, 1 ); 
			glEnable( GL_TEXTURE_2D ); 
			glBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA ); 
			glEnable( GL_BLEND ); 
			
			glTexImage2D( GL_TEXTURE_2D, 0, 4, GRIDSIZEY, GRIDSIZEX, 0, GL_RGBA, GL_UNSIGNED_BYTE, &SpatailGradientFlowX_dxTex[0][0][0] );  
					
			glBegin( GL_QUADS );     
				glTexCoord2f( 0., 0. ); 
				glVertex2f( 0., .5 ); 
				glTexCoord2f( 0., 1. ); 
				glVertex2f( 1., .5 ); 
				glTexCoord2f( 1., 1. ); 
				glVertex2f( 1., 1. ); 
				glTexCoord2f( 1., 0. ); 
				glVertex2f( 0., 1. ); 
			glEnd(); 

			glTexImage2D( GL_TEXTURE_2D, 0, 4, GRIDSIZEY, GRIDSIZEX, 0, GL_RGBA, GL_UNSIGNED_BYTE, &SpatailGradientFlowX_dyTex[0][0][0] );  					
		
			glBegin( GL_QUADS );     
				glTexCoord2f( 0., 0. ); 
				glVertex2f( 0., 0. ); 
				glTexCoord2f( 0., 1. ); 
				glVertex2f( 1., 0. ); 
				glTexCoord2f( 1., 1. ); 
				glVertex2f( 1., .5 ); 
				glTexCoord2f( 1., 0. ); 
				glVertex2f( 0., .5 ); 
			glEnd(); 

			glTexImage2D( GL_TEXTURE_2D, 0, 4, GRIDSIZEY, GRIDSIZEX, 0, GL_RGBA, GL_UNSIGNED_BYTE, &SpatailGradientFlowY_dxTex[0][0][0] );  
					
			glBegin( GL_QUADS );     
				glTexCoord2f( 0., 0. ); 
				glVertex2f( 1., .5 ); 
				glTexCoord2f( 0., 1. ); 
				glVertex2f( 2., .5 ); 
				glTexCoord2f( 1., 1. ); 
				glVertex2f( 2., 1. ); 
				glTexCoord2f( 1., 0. ); 
				glVertex2f( 1., 1. ); 
			glEnd(); 

			glTexImage2D( GL_TEXTURE_2D, 0, 4, GRIDSIZEY, GRIDSIZEX, 0, GL_RGBA, GL_UNSIGNED_BYTE, &SpatailGradientFlowY_dyTex[0][0][0] );  					
		
			glBegin( GL_QUADS );     
				glTexCoord2f( 0., 0. ); 
				glVertex2f( 1., 0. ); 
				glTexCoord2f( 0., 1. ); 
				glVertex2f( 2., 0. ); 
				glTexCoord2f( 1., 1. ); 
				glVertex2f( 2., .5 ); 
				glTexCoord2f( 1., 0. ); 
				glVertex2f( 1., .5 ); 
			glEnd(); 


			glDisable( GL_TEXTURE_2D ); 
				

		}

		if(FTLEFieldOn){
			glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP ); 
			glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP ); 
			glTexEnvf( GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE ); 
			filter = GL_LINEAR; 
		

			glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, filter ); 
			glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, filter ); 

			glPixelStorei( GL_UNPACK_ALIGNMENT, 1 ); 
			glEnable( GL_TEXTURE_2D ); 
			glBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA ); 
			//glEnable( GL_BLEND ); 
			
			glTexImage2D( GL_TEXTURE_2D, 0, 4, GRIDSIZEY, GRIDSIZEX, 0, GL_RGBA, GL_UNSIGNED_BYTE, &FTLEFieldTex[0][0][0] );  
					
			//glBegin( GL_QUADS );     
			//	glTexCoord2f( 0., 0. ); 
			//	glVertex2f( 2., 1. ); 
			//	glTexCoord2f( 0., 1. ); 
			//	glVertex2f( 0., 1. ); 
			//	glTexCoord2f( 1., 1. ); 
			//	glVertex2f( 0., 0. ); 
			//	glTexCoord2f( 1., 0. ); 
			//	glVertex2f( 2., 1. ); 
			//glEnd(); 

			glBegin( GL_QUADS );     
				glTexCoord2f( 0., 0. ); 
				glVertex2f( 0., 0. ); 
				glTexCoord2f( 0., 1. ); 
				glVertex2f( 2., 0. ); 
				glTexCoord2f( 1., 1. ); 
				glVertex2f( 2., 1. ); 
				glTexCoord2f( 1., 0. ); 
				glVertex2f( 0., 1. ); 
			glEnd(); 

			glDisable( GL_TEXTURE_2D ); 		

		}
	}

}


void Display_Model(void)
{
	GLint viewport[4];
	int jitter;

	glClearColor (1.0, 1.0, 1.0, 1.0);  // background for rendering color coding and lighting
	glGetIntegerv (GL_VIEWPORT, viewport);
 
	set_view(GL_RENDER);
	set_scene(GL_RENDER);
	display_shape(GL_RENDER);
	glFlush();
	glutSwapBuffers();
	glFinish();
}


void 
get_vector_field( float x, float y, float z, float t, float &vxp, float &vyp, float &vzp ) 
{ 
	//a(t) and b(t)
	float a_t = EPSILON * sin( OMEGA * PI * t);
	float b_t = 1.0f - 2.0f * EPSILON * sin( OMEGA * PI * t); 

	//f(x,t)
	float f_x = a_t * pow(x, 2) + b_t * x;

	//f'(x,t)
	float df_x = 2.0f * a_t * x + b_t;

	//u, v
	vxp = -1.0f * PI * AAAA * sin(PI * f_x) * cos ( PI * y);
	vyp = PI * AAAA * cos( PI * f_x) * sin( PI * y) * df_x;
	vzp = 0;

} 

void
get_vector_fieldF( float x, float y, float z, int t, float &vxp, float &vyp, float &vzp ) 
{ 
	/*printf("%f, %f\n", x, z);getchar();*/
	float dis = 0., tempdis = 0., tempvx = 0., tempzx = 0., tempvy = 0.;
	int flag = 0;
	for(int i= 0; i < vectorFieldQIL[t].size(); i++){
		dis = sqrt(pow(vectorFieldQIL[t][i]->tail->x/scaleqil - x,2) + pow(vectorFieldQIL[t][i]->tail->y/scaleqil - y,2) + pow(vectorFieldQIL[t][i]->tail->z/scaleqil - z,2));
		//printf("%f\n", dis);
		
		if(dis < 0.05){
			if(flag = 0){
				tempdis = dis;
				flag = 1;
			}
			vxp = vectorFieldQIL[t][i]->vect->x/scaleqil;
			vyp = 0.;
			vzp = vectorFieldQIL[t][i]->vect->z/scaleqil;
			
				if(dis < tempdis){
					vxp = vectorFieldQIL[t][i]->vect->x/scaleqil;
					vyp = 0.;
					vzp = vectorFieldQIL[t][i]->vect->z/scaleqil;
					tempdis = dis;

				/*	printf("%f\n", dis);	
					getchar();*/
				}
		}
	}

} 

void 
create_vector_field(float time){
	vector<node*> Ynodes;
	
	float incrementX, incrementY;
	incrementX = (float)2.0f/(float)GRIDSIZEX;
	incrementY = (float)1.0f/(float)GRIDSIZEY;
	
	vectorField.listOfNodes.clear();
	float initX, initY;

	initX = 0.0f;		
	for(int i = 0; i < GRIDSIZEX; i++){				
		initY = 0.0f;	
		Ynodes.clear();
		for(int j = 0; j < GRIDSIZEY; j++){					
			node* tempNode = (node*)malloc(sizeof(struct node));

			tempNode->x = initX;
			tempNode->y = initY;
			tempNode->z = 0.0f;
			/*printf("%d, %d, %d, %f, %f, %f\n", i, j, k, tempNode->x, tempNode->y, tempNode->z);getchar();*/
			get_vector_field(initX, initY, 0.0f, time, tempNode->vx, tempNode->vy, tempNode->vz);
			tempNode->magnitude = sqrt(pow(tempNode->vx ,2)+pow(tempNode->vy ,2)+pow(tempNode->vz ,2));
			Ynodes.push_back(tempNode);				
				
			initY += incrementY;
		}		
		vectorField.listOfNodes.push_back(Ynodes);
		initX += incrementX;
	}
	float max_magnitude = vectorField.listOfNodes[0][0]->magnitude, min_magnitude = vectorField.listOfNodes[0][0]->magnitude;
	for (int i = 0; i < GRIDSIZEX; i++){
		for (int j = 0; j < GRIDSIZEY; j++){			
			if (max_magnitude < vectorField.listOfNodes[i][j]->magnitude) max_magnitude = vectorField.listOfNodes[i][j]->magnitude;
			if (min_magnitude > vectorField.listOfNodes[i][j]->magnitude) min_magnitude = vectorField.listOfNodes[i][j]->magnitude;				
		}
	}

	float hsv[3], rgb[3];
	hsv[0] = 0.0f;rgb[0] = 0.0f;
	hsv[1] = 1.0f;rgb[1] = 0.0f;
	hsv[2] = 1.0f;rgb[2] = 0.0f;

	for(int i = 0; i < GRIDSIZEX; i++){
		for(int j = 0; j < GRIDSIZEY; j++){	
			hsv[0] = 240.0f - 240.0f*((vectorField.listOfNodes[i][j]->magnitude - min_magnitude)/(max_magnitude - min_magnitude));
			HsvRgb(hsv, rgb);
			vectorField.listOfNodes[i][j]->r = rgb[0];
			vectorField.listOfNodes[i][j]->g = rgb[1];
			vectorField.listOfNodes[i][j]->b = rgb[2];	
		}
	}		

}

void 
CalculateFlowMaps(){
	float incrementX, incrementY;
	incrementX = (float)2.0f/(float)GRIDSIZEX;
	incrementY = (float)1.0f/(float)GRIDSIZEY;
	
	float initX, initY;
	//Deploy Particles on Grid
	initX = 0.0f;		
	for(int i = 0; i < GRIDSIZEX; i++){				
		initY = 0.0f;			
		for(int j = 0; j < GRIDSIZEY; j++){		
			flowMapX[0][i][j] = initX;
			flowMapY[0][i][j] = initY;							
			initY += incrementY;
		}				
		initX += incrementX;
	}

	//Advect them using 4th order integration
	for(int t = 1; t < INTEGRATIONTIME*10; t++){
		for(int i = 0; i< GRIDSIZEX; i++){
			for(int j = 0; j< GRIDSIZEY; j++){
				//Initial Time
				float t1 = (t - 1.0f)/10.;

				//Initial Position
				float initialPositionX = flowMapX[(int)t1][i][j];
				float initialPositionY = flowMapY[(int)t1][i][j];
				float vectorMagnitude = 1.;
				//get the vector at t1
				float vectx = 0.0f, vecty = 0.0f, vectz = 0.0f;
				get_vector_field(initialPositionX, initialPositionY, 0.0f, t1, vectx, vecty, vectz);
				/*vectorMagnitude = sqrt(pow(vectx,2) + pow(vecty,2) + pow(vectz,2));*/
				
				float tempPositionX, tempPositionY;
				//Advect 
				if(vectorMagnitude !=0){
					tempPositionX = initialPositionX + vectx *(0.25f) * (1.0f/(vectorMagnitude/**(float)GRIDSIZEX*/));
					tempPositionY = initialPositionY + vecty *(0.25f) * (1.0f/(vectorMagnitude/**(float)GRIDSIZEY*/));
				}
				else{
					tempPositionX = initialPositionX + 0.0f;
					tempPositionY = initialPositionY + 0.0f;
				}			
				
				for (int k = 3; k > 0; k--){
					t1 = t - (0.25*(k))/10.;
					get_vector_field(tempPositionX, tempPositionY, 0.0f, t1, vectx, vecty, vectz);
					/*vectorMagnitude = sqrt(pow(vectx,2) + pow(vecty,2) + pow(vectz,2));*/

					if(vectorMagnitude !=0){
					tempPositionX = tempPositionX + vectx *(0.25f) * (1.0f/(vectorMagnitude/**(float)GRIDSIZEX*/));
					tempPositionY = tempPositionY + vecty *(0.25f) * (1.0f/(vectorMagnitude/**(float)GRIDSIZEY*/));
					}
					else{
						tempPositionX = tempPositionX + 0.0f;
						tempPositionY = tempPositionY + 0.0f;
					}		
					/*printf("%f, %f, %f\n", t1, tempPositionX, tempPositionY);*/					
				}
				//printf("%f, %f, %f", vectorMagnitude, tempPositionX, tempPositionY);getchar();
				/*printf("-------------------------\n");getchar();*/
				flowMapX[t][i][j] = tempPositionX;
				flowMapY[t][i][j] = tempPositionY;
				
			}
		}
	}
}

void 
CalculateFlowMapsRK4(){
	float incrementX, incrementY;
	incrementX = (float)2.0f/(float)GRIDSIZEX;
	incrementY = (float)1.0f/(float)GRIDSIZEY;
	
	float initX, initY;
	//Deploy Particles on Grid
	initX = 0.0f;		
	for(int i = 0; i < GRIDSIZEX; i++){				
		initY = 0.0f;			
		for(int j = 0; j < GRIDSIZEY; j++){		
			flowMapX[0][i][j] = initX;
			flowMapY[0][i][j] = initY;							
			initY += incrementY;
		}				
		initX += incrementX;
	}

	//Advect them using Runge Kutta 4th order integration
	for(int t = 1; t < INTEGRATIONTIME*10; t++){
		//getchar();
		for(int i = 0; i< GRIDSIZEX; i++){
			for(int j = 0; j< GRIDSIZEY; j++){
				//Initial Time
				float t1 = (t - 1.0f)/10.;
				//printf("%f, %d\n", t1, t);
				/*float scalingx = (float)GRIDSIZEX;
				float scalingy = (float)GRIDSIZEY;*/
				float scalingx = 80.;
				float scalingy = 80.;
				float vectorMagnitude = 1.;
				//Initial Position
				float initialPositionX = flowMapX[t - 1][i][j];
				float initialPositionY = flowMapY[t - 1][i][j];

				//get the vector at t1
				float vectx = 0.0f, vecty = 0.0f, vectz = 0.0f;
				get_vector_field(initialPositionX, initialPositionY, 0.0f, t1, vectx, vecty, vectz);
				vectorMagnitude = sqrt(pow(vectx,2) + pow(vecty,2) + pow(vectz,2));
				
				//vector a
				/*float ax = vectx/(vectorMagnitude*(float)GRIDSIZEX);
				float ay = vecty/(vectorMagnitude*(float)GRIDSIZEY);*/
				float ax = vectx/(vectorMagnitude*scalingx);
				float ay = vecty/(vectorMagnitude*scalingy);
				/*float ax = (vectx*vectorMagnitude)/scalingx;
				float ay = (vecty*vectorMagnitude)/scalingy;*/

				get_vector_field(initialPositionX + ax/2., initialPositionY + ay/2. , 0.0f, t1, vectx, vecty, vectz);
				vectorMagnitude = sqrt(pow(vectx,2) + pow(vecty,2) + pow(vectz,2));

				//vector b
				/*float bx = vectx/(vectorMagnitude*(float)GRIDSIZEX);
				float by = vecty/(vectorMagnitude*(float)GRIDSIZEY);*/
				float bx = vectx/(vectorMagnitude*scalingx);
				float by = vecty/(vectorMagnitude*scalingy);
				/*float bx = (vectx*vectorMagnitude)/scalingx;
				float by = (vecty*vectorMagnitude)/scalingy;*/

				get_vector_field(initialPositionX + bx/2., initialPositionY + by/2. , 0.0f, t1, vectx, vecty, vectz);
				vectorMagnitude = sqrt(pow(vectx,2) + pow(vecty,2) + pow(vectz,2));
				
				//vector c
				/*float cx = vectx/(vectorMagnitude*(float)GRIDSIZEX);
				float cy = vecty/(vectorMagnitude*(float)GRIDSIZEY);*/
				float cx = vectx/(vectorMagnitude*scalingx);
				float cy = vecty/(vectorMagnitude*scalingy);
				/*float cx = (vectx*vectorMagnitude)/scalingx;
				float cy = (vecty*vectorMagnitude)/scalingy;*/

				get_vector_field(initialPositionX + cx, initialPositionY + cy , 0.0f, t1, vectx, vecty, vectz);
				vectorMagnitude = sqrt(pow(vectx,2) + pow(vecty,2) + pow(vectz,2));

				//vector d
				/*float dx = vectx/(vectorMagnitude*(float)GRIDSIZEX);
				float dy = vecty/(vectorMagnitude*(float)GRIDSIZEY);*/
				float dx = vectx/(vectorMagnitude*scalingx);
				float dy = vecty/(vectorMagnitude*scalingy);
				/*float dx = (vectx*vectorMagnitude)/scalingx;
				float dy = (vecty*vectorMagnitude)/scalingy;*/

				float tempx= initialPositionX + (ax + 2.0f*bx + 2.0f*cx + dx)/6.0f ;
				float tempy= initialPositionY + (ay + 2.0f*by + 2.0f*cy + dy)/6.0f ;

				flowMapX[t][i][j] = tempx;
				flowMapY[t][i][j] = tempy;
				
				/*if(tempx >= 0. && tempx <=2. && tempy >= 0. && tempy <= 1.){
					flowMapX[t][i][j] = tempx;
					flowMapY[t][i][j] = tempy;
				}
				else{
					flowMapX[t][i][j] = initialPositionX;
					flowMapY[t][i][j] = initialPositionY;
				}*/

			}
		}
	}
}


void 
CalculateSpatialGradients(){
	for (int t = 1; t < INTEGRATIONTIME*10; t++){
		for (int i = 1; i < GRIDSIZEX - 1; i++){
			for (int j = 1; j < GRIDSIZEY - 1; j++){

				SpatailGradientFlowX_dx[t][i][j] = (flowMapX[t][i+1][j] - flowMapX[t][i-1][j])/(flowMapX[0][i+1][j] - flowMapX[0][i-1][j]);
				SpatailGradientFlowX_dy[t][i][j] = (flowMapX[t][i][j+1] - flowMapX[t][i][j-1])/(flowMapY[0][i][j+1] - flowMapY[0][i][j-1]);
				SpatailGradientFlowY_dx[t][i][j] = (flowMapY[t][i+1][j] - flowMapY[t][i-1][j])/(flowMapX[0][i+1][j] - flowMapX[0][i-1][j]);
				SpatailGradientFlowY_dy[t][i][j] = (flowMapY[t][i][j+1] - flowMapY[t][i][j-1])/(flowMapY[0][i][j+1] - flowMapY[0][i][j-1]);							
				//printf("%f, %f\n", SpatailGradientFlowX_dy[t][i][j], SpatailGradientFlowY_dx[t][i][j]);
				
				/*SpatailGradientFlowX_dx[t][i][j] = (flowMapX[t][i+1][j] - flowMapX[t][i-1][j])/(flowMapX[t-1][i+1][j] - flowMapX[t-1][i-1][j]);
				SpatailGradientFlowX_dy[t][i][j] = (flowMapX[t][i][j+1] - flowMapX[t][i][j-1])/(flowMapY[t-1][i][j+1] - flowMapY[t-1][i][j-1]);
				SpatailGradientFlowY_dx[t][i][j] = (flowMapY[t][i+1][j] - flowMapY[t][i-1][j])/(flowMapX[t-1][i+1][j] - flowMapX[t-1][i-1][j]);
				SpatailGradientFlowY_dy[t][i][j] = (flowMapY[t][i][j+1] - flowMapY[t][i][j-1])/(flowMapY[t-1][i][j+1] - flowMapY[t-1][i][j-1]);*/
			}
		}/*if(t>2)getchar();*/
	}
}
void 
CalculateFTLEField(){

	float Phi[2][2];float** CachyGreenDistrotionTensor;
	float eigenVector;
	for (int t = 1; t < INTEGRATIONTIME*10; t++){
		for (int i = 1; i < GRIDSIZEX - 1; i++){
			for (int j = 1; j < GRIDSIZEY - 1; j++){
				Phi[0][0] = SpatailGradientFlowX_dy[t][i][j];								
				Phi[0][1] = SpatailGradientFlowX_dx[t][i][j];								
				Phi[1][0] = SpatailGradientFlowY_dy[t][i][j];								
				Phi[1][1] = SpatailGradientFlowY_dx[t][i][j];		
				CachyGreenDistrotionTensor = CalculateCauchyGreen(Phi);
				/*printf("%d, %d, %f, %f, %f, %f\n", i, j, CachyGreenDistrotionTensor[0][0], CachyGreenDistrotionTensor[0][1], CachyGreenDistrotionTensor[1][0], CachyGreenDistrotionTensor[1][1]); */
				eigenVector = CalcualteEigenVector(CachyGreenDistrotionTensor);
				
				if(eigenVector > 0.){
					FTLEfield[t][i][j] = 10.0f * log(eigenVector)/(float)t;
				}
				else{
					/*FTLEfield[t][i][j] = FTLEfield[t-1][i][j];*/
					FTLEfield[t][i][j] = 0.;
				}
				//printf("%d, %d, %f, %f\n", i, j, eigenVector, FTLEfield[t][i][j]);
			}
		}
		//getchar();
	}


}

float 
CalcualteEigenVector(float** inputMatrix){
	float a = 1.;
	float b = -(inputMatrix[0][0] + inputMatrix[1][1]);
	float c = inputMatrix[0][0]*inputMatrix[1][1] - inputMatrix[0][1]*inputMatrix[1][0];

	float det = pow(b,2) - 4*a*c;
	float lambda1, lambda2;
	if(det >= 0){
		lambda1 = (-1.*b + sqrt(det))/(2.*a);
		lambda2 = (-1.*b - sqrt(det))/(2.*a);
	}
	else{
		/*printf("det < 0 ----%f\n", det);*/
		return 0.0f;
	}
	//printf("%f,%f\n",lambda1,lambda2);
	if (lambda1>lambda2) return lambda1;
	
	return lambda2;

}

float** 
CalculateCauchyGreen(float inputMatrix[][2]){
	float** CauchyGreenDistortion;
	CauchyGreenDistortion = new float*[2];	
	CauchyGreenDistortion[0] = new float[2];
	CauchyGreenDistortion[1] = new float[2];
	//Phi' * Phi
	CauchyGreenDistortion[0][0] = inputMatrix[0][0]*inputMatrix[0][0]+inputMatrix[1][0]*inputMatrix[1][0];
	CauchyGreenDistortion[0][1] = inputMatrix[0][0]*inputMatrix[0][1]+inputMatrix[1][0]*inputMatrix[1][1];
	CauchyGreenDistortion[1][0] = inputMatrix[0][1]*inputMatrix[0][0]+inputMatrix[1][1]*inputMatrix[1][0];
	CauchyGreenDistortion[1][1] = inputMatrix[0][1]*inputMatrix[0][1]+inputMatrix[1][1]*inputMatrix[1][1];
	/*CauchyGreenDistortion[0][0] = inputMatrix[0][0]*inputMatrix[0][0]+inputMatrix[0][1]*inputMatrix[0][1];
	CauchyGreenDistortion[0][1] = inputMatrix[0][0]*inputMatrix[1][0]+inputMatrix[0][1]*inputMatrix[1][1];
	CauchyGreenDistortion[1][0] = inputMatrix[1][0]*inputMatrix[0][0]+inputMatrix[1][1]*inputMatrix[0][1];
	CauchyGreenDistortion[1][1] = inputMatrix[1][0]*inputMatrix[1][0]+inputMatrix[1][1]*inputMatrix[1][1];*/

			//printf("%f,",inputMatrix[i][j]);

	return CauchyGreenDistortion;

}



void 
genSpatialGradientTexture(int inputTime)
{

	float maxxx, minxx, maxyx, minyx, maxxy, minxy, maxyy, minyy;
	float hsv[3], rgb[3];

	if(inputTime < INTEGRATIONTIME*10){
		maxxx = SpatailGradientFlowX_dx[inputTime][0][0]; minxx = SpatailGradientFlowX_dx[inputTime][0][0];
		maxxy = SpatailGradientFlowX_dy[inputTime][0][0]; minxx = SpatailGradientFlowX_dy[inputTime][0][0];
		maxyx = SpatailGradientFlowY_dx[inputTime][0][0]; minyx = SpatailGradientFlowY_dx[inputTime][0][0];
		maxyy = SpatailGradientFlowY_dy[inputTime][0][0]; minyy = SpatailGradientFlowY_dy[inputTime][0][0];
		
		for(int i = 0; i < GRIDSIZEX; i++){
			for (int j = 0; j < GRIDSIZEY; j++){
				if (maxxx < SpatailGradientFlowX_dx[inputTime][i][j]) maxxx = SpatailGradientFlowX_dx[inputTime][i][j];
				if (minxx > SpatailGradientFlowX_dx[inputTime][i][j]) minxx = SpatailGradientFlowX_dx[inputTime][i][j];	

				if (maxxy < SpatailGradientFlowX_dy[inputTime][i][j]) maxxy = SpatailGradientFlowX_dy[inputTime][i][j];
				if (minxy > SpatailGradientFlowX_dy[inputTime][i][j]) minxy = SpatailGradientFlowX_dy[inputTime][i][j];

				if (maxyx < SpatailGradientFlowY_dx[inputTime][i][j]) maxyx = SpatailGradientFlowY_dx[inputTime][i][j];
				if (minyx > SpatailGradientFlowY_dx[inputTime][i][j]) minyx = SpatailGradientFlowY_dx[inputTime][i][j];

				if (maxyy < SpatailGradientFlowY_dy[inputTime][i][j]) maxyy = SpatailGradientFlowY_dy[inputTime][i][j];
				if (minyy > SpatailGradientFlowY_dy[inputTime][i][j]) minyy = SpatailGradientFlowY_dy[inputTime][i][j];

			}
		}



		for(int i = 0; i < GRIDSIZEX; i++){
			for (int j = 0; j < GRIDSIZEY; j++){				
					/*TextureXY[zz][y][x][0] = (unsigned char) ( 255.*r + .5 ); 
					TextureXY[zz][y][x][1] = (unsigned char) ( 255.*g + .5 ); 
					TextureXY[zz][y][x][2] = (unsigned char) ( 255.*b + .5 ); 
					TextureXY[zz][y][x][3] = (unsigned char) ( 255.*alpha + .5 ); */
					hsv[0] = 0.0f;rgb[0] = 0.0f;
					hsv[1] = 1.0f;rgb[1] = 0.0f;
					hsv[2] = 1.0f;rgb[2] = 0.0f;

					hsv[0] = 240.0f - 240.0f*((SpatailGradientFlowX_dx[inputTime][i][j] - minxx)/(maxxx - minxx));
					HsvRgb(hsv, rgb);

					SpatailGradientFlowX_dxTex[i][j][0] = (unsigned char) ( 255.*rgb[0] + .5 ); 					
					SpatailGradientFlowX_dxTex[i][j][1] = (unsigned char) ( 255.*rgb[1] + .5 ); 					
					SpatailGradientFlowX_dxTex[i][j][2] = (unsigned char) ( 255.*rgb[2] + .5 ); 					
					/*SpatailGradientFlowX_dxTex[i][j][3] = (unsigned char) ( 255.*ALPHA + .5 ); 		*/
					SpatailGradientFlowX_dxTex[i][j][3] = 255.; 		

					hsv[0] = 0.0f;rgb[0] = 0.0f;
					hsv[1] = 1.0f;rgb[1] = 0.0f;
					hsv[2] = 1.0f;rgb[2] = 0.0f;

					hsv[0] = 240.0f - 240.0f*((SpatailGradientFlowX_dy[inputTime][i][j] - minxy)/(maxxy - minxy));
					HsvRgb(hsv, rgb);

					SpatailGradientFlowX_dyTex[i][j][0] = (unsigned char) ( 255.*rgb[0] + .5 ); 					
					SpatailGradientFlowX_dyTex[i][j][1] = (unsigned char) ( 255.*rgb[1] + .5 ); 					
					SpatailGradientFlowX_dyTex[i][j][2] = (unsigned char) ( 255.*rgb[2] + .5 ); 					
					/*SpatailGradientFowX_dyTex[i][j][3] = (unsigned char) ( 255.*ALPHA + .5 ); 	*/
					SpatailGradientFlowX_dyTex[i][j][3] = 255.; 	

					hsv[0] = 0.0f;rgb[0] = 0.0f;
					hsv[1] = 1.0f;rgb[1] = 0.0f;
					hsv[2] = 1.0f;rgb[2] = 0.0f;

					hsv[0] = 240.0f - 240.0f*((SpatailGradientFlowY_dx[inputTime][i][j] - minyx)/(maxyx - minyx));
					HsvRgb(hsv, rgb);

					SpatailGradientFlowY_dxTex[i][j][0] = (unsigned char) ( 255.*rgb[0] + .5 ); 					
					SpatailGradientFlowY_dxTex[i][j][1] = (unsigned char) ( 255.*rgb[1] + .5 ); 					
					SpatailGradientFlowY_dxTex[i][j][2] = (unsigned char) ( 255.*rgb[2] + .5 ); 					
					//SpatailGradientFlowY_dxTex[i][j][3] = (unsigned char) ( 255.*ALPHA + .5 ); 	
					SpatailGradientFlowY_dxTex[i][j][3] = 255.; 	

					hsv[0] = 0.0f;rgb[0] = 0.0f;
					hsv[1] = 1.0f;rgb[1] = 0.0f;
					hsv[2] = 1.0f;rgb[2] = 0.0f;

					hsv[0] = 240.0f - 240.0f*((SpatailGradientFlowY_dy[inputTime][i][j] - minyy)/(maxyy - minyy));
					HsvRgb(hsv, rgb);

					SpatailGradientFlowY_dyTex[i][j][0] = (unsigned char) ( 255.*rgb[0] + .5 ); 					
					SpatailGradientFlowY_dyTex[i][j][1] = (unsigned char) ( 255.*rgb[1] + .5 ); 					
					SpatailGradientFlowY_dyTex[i][j][2] = (unsigned char) ( 255.*rgb[2] + .5 ); 					
					/*SpatailGradientFlowY_dyTex[i][j][3] = (unsigned char) ( 255.*ALPHA + .5 ); 	*/
					SpatailGradientFlowY_dyTex[i][j][3] = 255.; 	
					
					
				
			}
		}
	}

}

void 
genFTLETexture(int inputTime)
{

	float max, min;
	float hsv[3], rgb[3];
	
	if(inputTime < INTEGRATIONTIME*10){
		max = FTLEfield[inputTime][0][0]; min = FTLEfield[inputTime][0][0];
		
		
		for(int i = 0; i < GRIDSIZEX; i++){
			for (int j = 0; j < GRIDSIZEY; j++){
				if (max < FTLEfield[inputTime][i][j]) max = FTLEfield[inputTime][i][j];
				if (min > FTLEfield[inputTime][i][j]) min = FTLEfield[inputTime][i][j];	
				

			}
			//printf("%f, %f", min, max);
		}



		for(int i = 0; i < GRIDSIZEX; i++){
			for (int j = 0; j < GRIDSIZEY; j++){				
					/*TextureXY[zz][y][x][0] = (unsigned char) ( 255.*r + .5 ); 
					TextureXY[zz][y][x][1] = (unsigned char) ( 255.*g + .5 ); 
					TextureXY[zz][y][x][2] = (unsigned char) ( 255.*b + .5 ); 
					TextureXY[zz][y][x][3] = (unsigned char) ( 255.*alpha + .5 ); */
					hsv[0] = 0.0f;rgb[0] = 0.0f;
					hsv[1] = 1.0f;rgb[1] = 0.0f;
					hsv[2] = 1.0f;rgb[2] = 0.0f;

					hsv[0] = 240.0f - 240.0f*((FTLEfield[inputTime][i][j] - min)/(max - min));
					HsvRgb(hsv, rgb);

					FTLEFieldTex[i][j][0] = (unsigned char) ( 255.*rgb[0] + .5 ); 					
					FTLEFieldTex[i][j][1] = (unsigned char) ( 255.*rgb[1] + .5 ); 					
					FTLEFieldTex[i][j][2] = (unsigned char) ( 255.*rgb[2] + .5 ); 					
					/*SpatailGradientFlowX_dxTex[i][j][3] = (unsigned char) ( 255.*ALPHA + .5 ); 		*/
					FTLEFieldTex[i][j][3] = (unsigned char) (255.); 							
				
			}
		}
	}

}



qilTriangelMesh* loaddata(){
	/*FILE * pFile;
	pFile = fopen("deteiledfloor.qil", "r");
	
	char mystring[100];
	puts(mystring);

	printf("%02X\n", mystring[0]);
*/
	vector<qilTriangle*> meshToReturn;
	vector<qilVertex*> vertexToreturn;
	qilTriangelMesh* triangleMeshToReturn= (qilTriangelMesh*) malloc(sizeof(qilTriangelMesh));

	FILE *this_file = fopen("dg48.ply", "r");
	poly = new Polyhedron (this_file);
	fclose(this_file);
	poly->initialize();
	for(int i = 0; i < poly->nverts; i++){
		/*printf("%f, %f, %f\n", poly->vlist[i]->x, poly->vlist[i]->y, poly->vlist[i]->z);*/
		qilVertex *tempVert1 = (qilVertex*)malloc(sizeof(struct qilVertex));
		tempVert1->x = poly->vlist[i]->x;
		tempVert1->y = poly->vlist[i]->y;
		tempVert1->z = poly->vlist[i]->z;
		vertexToreturn.push_back(tempVert1);

	}
	for(int i = 0; i < poly->ntris ; i++){
		/*printf("%d, %d, %d\n", poly->tlist[i]->verts[0]->index, poly->tlist[i]->verts[1]->index, poly->tlist[i]->verts[2]->index);*/
		qilTriangle* temp = (qilTriangle*)malloc(sizeof(qilTriangle*));
		temp->verts[0] = poly->tlist[i]->verts[0]->index;
		temp->verts[1] = poly->tlist[i]->verts[1]->index;
		temp->verts[2] = poly->tlist[i]->verts[2]->index;
		meshToReturn.push_back(temp);
	}
	//triangleMeshToReturn->vertices = (qilVertex*) malloc(sizeof(struct qilVertex));
	triangleMeshToReturn->vertices = vertexToreturn;
	triangleMeshToReturn->triangles = meshToReturn;

	/*for(int i = 0; i < triangleMeshToReturn->triangles.size() ; i++){
		triangleMeshToReturn->triangles[i]->centroid->x = (triangleMeshToReturn->vertices[triangleMeshToReturn->triangles[i]->verts[0]]->x + triangleMeshToReturn->vertices[triangleMeshToReturn->triangles[i]->verts[1]]->x + triangleMeshToReturn->vertices[triangleMeshToReturn->triangles[i]->verts[2]]->x)/3;
		triangleMeshToReturn->triangles[i]->centroid->y = (triangleMeshToReturn->vertices[triangleMeshToReturn->triangles[i]->verts[0]]->y + triangleMeshToReturn->vertices[triangleMeshToReturn->triangles[i]->verts[1]]->y + triangleMeshToReturn->vertices[triangleMeshToReturn->triangles[i]->verts[2]]->y)/3;
		triangleMeshToReturn->triangles[i]->centroid->z = (triangleMeshToReturn->vertices[triangleMeshToReturn->triangles[i]->verts[0]]->z + triangleMeshToReturn->vertices[triangleMeshToReturn->triangles[i]->verts[1]]->z + triangleMeshToReturn->vertices[triangleMeshToReturn->triangles[i]->verts[2]]->z)/3;
	}*/

	
	return triangleMeshToReturn;
}

void calculateCentroids(){
	for(int i = 0; i < detailedFloor->triangles.size() ; i++){
		detailedFloor->triangles[i]->centroid = (qilVertex*)malloc(sizeof(struct qilVertex));
		detailedFloor->triangles[i]->centroid->x = (detailedFloor->vertices[detailedFloor->triangles[i]->verts[0]]->x + detailedFloor->vertices[detailedFloor->triangles[i]->verts[1]]->x + detailedFloor->vertices[detailedFloor->triangles[i]->verts[2]]->x)/3;
		detailedFloor->triangles[i]->centroid->y = (detailedFloor->vertices[detailedFloor->triangles[i]->verts[0]]->y + detailedFloor->vertices[detailedFloor->triangles[i]->verts[1]]->y + detailedFloor->vertices[detailedFloor->triangles[i]->verts[2]]->y)/3;
		detailedFloor->triangles[i]->centroid->z = (detailedFloor->vertices[detailedFloor->triangles[i]->verts[0]]->z + detailedFloor->vertices[detailedFloor->triangles[i]->verts[1]]->z + detailedFloor->vertices[detailedFloor->triangles[i]->verts[2]]->z)/3;
	}

}

void loadVector(){	
	/*
	float a1[200] = { 418,1312.4,38,1795.78,1318.53,38,1787.34,1272.8,38,1789.11,1280.95,38,1781.65,1339.02,38,1750.69,1344.21,38,1744,1408.37,38,1689.42,1412.18,38,1684.23,1444.9,38,1671.37,1451.83,38,1663.08,1377.81,38,1678.19,1367.92,38,1688.52,1355.84,38,1674.21,1347.59,38,1682.04,1327.13,38,1673.15,1322.75,38,1678.19,1306.53,38,1670.16,1301.14,38,1676.45,1277.84,38,1669.57,1269.56,38,1678.24,1266.97,38,1720.15,1258.02,38,1732.92,1215.42,38,1785.54,1209.02,38,1793.74,1184.92,38,1790.58,1176.04,38,1801.17,1168.34,38,1786.03,1161.67,38,1795.15};
	float a2[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a3[200] = { 424,1277.9,38,1825.47,1286.49,38,1816.53,1266.8,38,1818.7,1269.77,38,1815.46,1312.36,38,1791.97,1317.77,38,1787.34,1273.82,38,1785.35,1279.3,38,1779.68,1340.6,38,1746.12,1349.88,38,1736.58,1412.7,38,1683.72,1417.57,38,1676.97,1379.68,38,1679.75,1370.41,38,1690.74,1356.99,38,1674.76,1349.04,38,1683.39,1321.76,38,1675.58,1315.54,38,1685.77,1302.32,38,1672.12,1296.54,38,1677.96,1283.97,38,1669.92,1277.76,38,1674.95,1266.43,38,1721.11,1260.99,38,1728.92,1211.19,38,1788.96,1199.68,38,1806.18,1181.88,38,1793.59,1179.7,38,1797.56,1168.73,38,1786.03,1165.08,38,1792.91 };
	float a4[200] = { 427,1314.28,38,1791.97,1318.53,38,1787.34,1283.23,38,1782.03,1285.29,38,1778.38,1342.8,38,1743.35,1348.86,38,1736.58,1415.52,38,1678.68,1425.56,38,1669.51,1379.52,38,1683.41,1372.28,38,1691.86,1356.12,38,1676.47,1346.53,38,1685.77,1325.57,38,1673.6,1319.19,38,1680.22,1300.43,38,1670.51,1295.91,38,1678.47,1280.51,38,1671.89,1275.9,38,1676.5,1259.07,38,1724.92,1247.74,38,1736.98,1217.83,38,1783.89,1208.12,38,1797.17,1183.92,38,1792.68,1177.18,38,1802.52,1203.87,38,1799.01,1199.85,38,1802.63,1168.63,38,1789.9,1164.18,38,1794.23 };
	float a5[200] = { 433,1284.74,38,1820.29,1293.41,38,1810.45,1268.12,38,1820.97,1272.77,38,1814.91,1313.14,38,1791.19,1323.15,38,1780.85,1282.88,38,1776.06,1284.93,38,1772.49,1342.46,38,1743.35,1347.84,38,1736.24,1422.41,38,1669.22,1431.09,38,1664.42,1371.9,38,1689.13,1361.29,38,1701.46,1349.82,38,1677.4,1345.45,38,1681.22,1329.29,38,1676.99,1325.51,38,1681.92,1302.14,38,1677.79,1297.64,38,1683.31,1280.09,38,1670.62,1275.76,38,1675.73,1261.67,38,1723.95,1250.52,38,1734.55,1218.91,38,1785.93,1207.33,38,1798.49,1181.08,38,1793.59,1175.95,38,1800.24,1163.06,38,1794.74,1160.82,38,1800.92 };
	float a6[200] = { 436,1281.69,38,1821.3,1288.33,38,1814.7,1263.2,38,1817.44,1267.76,38,1812.72,1312.36,38,1791.97,1321.56,38,1781.32,1280.2,38,1778.83,1291.01,38,1769.72,1345.16,38,1745.77,1351.22,38,1740.59,1422.96,38,1672.59,1427.6,38,1666.43,1371.26,38,1688.33,1367.14,38,1694.23,1351.3,38,1678.46,1343.14,38,1685.95,1330.73,38,1677.55,1324.46,38,1683.5,1304.51,38,1673.94,1298.01,38,1680.92,1285.06,38,1670.47,1275.71,38,1679.16,1260.62,38,1725.61,1252.9,38,1736.65,1208.9,38,1795.85,1199.23,38,1810.73,1183.02,38,1795.3,1177.65,38,1801.17,1162.95,38,1796.07,1156.93,38,1803.65,1196.87,38,1804.41,1188.95,38,1812.71 };
	float a7[200] = { 439,1283.03,38,1824.45,1293.47,38,1813.18,1266.55,38,1812.72,1271.72,38,1805.91,1316.36,38,1785.35,1323.98,38,1779.28,1286.26,38,1776.43,1291.43,38,1772.06,1354.32,38,1740.93,1362.16,38,1731.31,1370.06,38,1688.82,1362.53,38,1698.92,1346,38,1683.09,1335.9,38,1695.16,1323.15,38,1680.83,1318.77,38,1688.54,1298.25,38,1679.04,1294.03,38,1686.22,1278.34,38,1678.68,1272.13,38,1686.11,1258.79,38,1724.3,1255.82,38,1731.22,1212.39,38,1790.64,1204.27,38,1802.09,1198.06,38,1806.18,1189.04,38,1819.68,1178.13,38,1799.82,1174.1,38,1809.7,1169.06,38,1795.54,1159.97,38,1803.22 };
	float a8[200] = { 442,1283.85,38,1824.45,1297.08,38,1806.4,1266.07,38,1808.21,1270.02,38,1803.36,1316.42,38,1789.83,1322.68,38,1781.69,1288.65,38,1773.24,1293.52,38,1768.12,1351.19,38,1739.95,1357.06,38,1734.27,1422.53,38,1664.46,1425.97,38,1660.13,1371.87,38,1688.33,1364.13,38,1701.78,1347.26,38,1682.82,1341.47,38,1689.17,1330.15,38,1681.74,1320.01,38,1692.68,1299.1,38,1682.24,1289.63,38,1694.21,1273.39,38,1679.64,1266.36,38,1690.2,1255.25,38,1729.22,1244.85,38,1740.1,1210.07,38,1794.54,1199.09,38,1811.58,1179.23,38,1799.3,1171.4,38,1808.45,1167.86,38,1795.54,1159.37,38,1804.05,1189.19,38,1808.66,1180.68,38,1819.87 };
	float a9[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a10[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
*/
	/*float a11[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a12[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a13[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a14[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a15[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a16[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a17[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a18[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a19[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a20[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };

	float a21[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a22[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a23[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a24[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a25[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a26[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a27[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a28[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a29[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a30[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };

	float a31[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a32[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a33[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a34[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a35[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a36[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a37[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a38[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a39[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
	float a40[200] = { 421,1312.9,38,1794.51,1319.95,38,1788.97,1279.96,38,1783.63,1286.78,38,1775.98,1342.71,38,1747.13,1346.64,38,1742.67,1375.9,38,1681.76,1368.86,38,1689.63,1412.6,38,1682.95,1422.07,38,1674.59,1351.05,38,1675.88,1340,38,1687.25,1326.6,38,1672.6,1322.61,38,1680.26,1306.73,38,1670.91,1299.59,38,1678.57,1281.9,38,1667.62,1276.54,38,1674.95,1269.82,38,1719.81,1262.87,38,1728.56,1213.11,38,1786.42,1207.69,38,1792.84,1187.97,38,1790.08,1185.16,38,1794.89,1170.39,38,1787.31,1168.54,38,1791.2 };
*/
	vectorFieldQIL.empty();
   char str[2000];
   FILE *fp;
 
   fp = fopen("vector1.csv", "r");
   
   while(fgets(str,sizeof(str),fp) != NULL)
   {
      // strip trailing '\n' if it exists
	  vector< vectorsamp* > vectorofStruct;
      int len = strlen(str)-1;
      if(str[len] == '\n') 
         str[len] = 0;
	  //printf("\n%d, %s\n\n", len, str);
	  char * pch;
	  pch = strtok (str,",");
	  float temp[6];

	  while (pch != NULL)
	  {
		//printf ("%s\n",pch);
		vectorsamp* tempVectSamp = (vectorsamp*) malloc(sizeof(struct vectorsamp));
		for(int i = 0; i < 6; i++){
		pch = strtok (NULL, ",");
		if(pch == NULL) break;
		temp[i] = atof(pch);
		}
		qilVertex* temp1 = (qilVertex*) malloc(sizeof(struct qilVertex));
		qilVertex* temp2 = (qilVertex*) malloc(sizeof(struct qilVertex));
		qilVertex* temp3 = (qilVertex*) malloc(sizeof(struct qilVertex));
		
		temp1->x = temp[0];
		temp1->y = 0.;
		temp1->z = temp[2];

		temp2->x = temp[3];
		temp2->y = 0.;
		temp2->z = temp[5];

		temp3->x = temp[3] - temp[0];
		temp3->y = 0.;
		temp3->z = temp[5] - temp[2];

		tempVectSamp->tail = temp1;
		tempVectSamp->head = temp2;
		tempVectSamp->vect = temp3;
		vectorofStruct.push_back(tempVectSamp);

	  }
	  vectorFieldQIL.push_back(vectorofStruct);
   }
   fclose(fp);


}

void createVectorFieldF(int time){	
	vector<node*> Ynodes;
	
	float incrementX, incrementY;

	incrementX = (float)(1514. - 936.)/((float)GRIDSIZEY*scaleqil);
	incrementY = (float)(2043. - 1147.)/((float)GRIDSIZEX*scaleqil);
	
	vectorField.listOfNodes.clear();
	float initX, initY;

	initX = 936.0f/scaleqil;		
	for(int i = 0; i < GRIDSIZEY; i++){		
		
		initY = 1147.0f/scaleqil;	
		Ynodes.clear();
		for(int j = 0; j < GRIDSIZEX; j++){					
			node* tempNode = (node*)malloc(sizeof(struct node));

			tempNode->x = initX;
			tempNode->y = 0.0;
			tempNode->z = initY;
			//printf("%d, %d, %f, %f, %f\n", i, j, tempNode->x, tempNode->y, tempNode->z);
			get_vector_fieldF(initX, 0.0f, initY, time, tempNode->vx, tempNode->vy, tempNode->vz);
			tempNode->magnitude = sqrt(pow(tempNode->vx ,2)+pow(tempNode->vy ,2)+pow(tempNode->vz ,2));
			Ynodes.push_back(tempNode);				
				
			initY += incrementY;
		}		
		/*if(i%10 == 0)getchar();*/
		vectorField.listOfNodes.push_back(Ynodes);
		initX += incrementX;
	}
	/*float max_magnitude = vectorField.listOfNodes[0][0]->magnitude, min_magnitude = vectorField.listOfNodes[0][0]->magnitude;
	for (int i = 0; i < GRIDSIZEX; i++){
		for (int j = 0; j < GRIDSIZEY; j++){			
			if (max_magnitude < vectorField.listOfNodes[i][j]->magnitude) max_magnitude = vectorField.listOfNodes[i][j]->magnitude;
			if (min_magnitude > vectorField.listOfNodes[i][j]->magnitude) min_magnitude = vectorField.listOfNodes[i][j]->magnitude;				
		}
	}

	float hsv[3], rgb[3];
	hsv[0] = 0.0f;rgb[0] = 0.0f;
	hsv[1] = 1.0f;rgb[1] = 0.0f;
	hsv[2] = 1.0f;rgb[2] = 0.0f;

	for(int i = 0; i < GRIDSIZEX; i++){
		for(int j = 0; j < GRIDSIZEY; j++){	
			hsv[0] = 240.0f - 240.0f*((vectorField.listOfNodes[i][j]->magnitude - min_magnitude)/(max_magnitude - min_magnitude));
			HsvRgb(hsv, rgb);
			vectorField.listOfNodes[i][j]->r = rgb[0];
			vectorField.listOfNodes[i][j]->g = rgb[1];
			vectorField.listOfNodes[i][j]->b = rgb[2];	
		}
	}		*/

}